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Abstract 

Consider a system governed by the time-dependent Schrodinger equation in its ground state. 
When subjected to weak (size e) parametric forcing by an "ionizing field" (time- varying), the 
state decays with advancing time due to coupling of the bound state to radiation modes. The 
decay-rate of this metastable state is governed by Fermi's Golden Rule, T[V], which depends on 
the potential V and the details of the forcing. We pose the potential design problem: find V opt 
which minimizes T[V] (maximizes the lifetime of the state) over an admissible class of potentials 
with fixed spatial support. We formulate this problem as a constrained optimization problem 
and prove that an admissible optimal solution exists. Then, using quasi-Newton methods, we 
compute locally optimal potentials. These have the structure of a truncated periodic poten- 
tial with a localized defect. In contrast to optimal structures for other spectral optimization 
problems, our optimizing potentials appear to be interior points of the constraint set and to be 
smooth. The multi-scale structures that emerge incorporate the physical mechanisms of energy 
confinement via material contrast and interference effects. 

An analysis of locally optimal potentials reveals local optimality is attained via two mech- 
anisms: (i) decreasing the density of states near a resonant frequency in the continuum and 
(ii) tuning the oscillations of extended states to make T[V], an oscillatory integral, small. Our 
approach achieves lifetimes, ~ (e 2 r[F]) \ for locally optimal potentials with <~ 0(1O 9 ) as 
compared with T -1 <~ O(10 2 ) for a typical potential. Finally, we explore the performance of 
optimal potentials via simulations of the time-evolution. 

1 Introduction 

In many problems of fundamental and applied science a scenario arises where, due to physical laws 
or engineering design, a state of the system is metastable; the state is long-lived but of finite lifetime 
due to coupling or leakage to an environment. In settings as diverse as linear and nonlinear optics, 
cavity QED, Bose-Einstein condensation (BEC), and quantum computation, one is interested in 
the manipulation of the lifetime of such metastable states. Our goal in this paper is to explore the 
problem of maximizing the lifetime of a metastable state for a class of ionization problems. The 
approach we take is applicable to a wide variety of linear and nonlinear problems. Specific examples 
where metastable states exist include: 
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1. the excited state of an atom, e.g. hydrogen, where due to coupling to a photon field, the 
atom in its excited state spontaneously undergoes a transition to its ground state, after some 
excited state lifetime [4], 

2. an approximate bound state (quasi-mode) of a quantum system, e.g. atom in a cavity or 
BEC in a magnetic trap, which leaks ("tunnels") out of the cavity and whose wave function 
decays with advancing time [28], 

3. an approximate guided mode of an optical waveguide which, due to scattering, bends in the 
waveguide, or diffraction, leaks out of the structure, resulting in attenuation of the wave-field 
within the waveguide with increasing propagation distance [24], and 

4. a "scatterer" which confines "rays" , but leaks energy to spatial infinity due to their wave 
nature, e.g. Helmholtz resonator, traps rays between obstacles [23]. 

These examples are representative of a class of extended (infinite spatial domain), yet energy- 
preserving (closed) systems, where the mechanism for energy loss is scattering loss, the escape of 
energy from a compact spatial domain to spatial infinity. 

Such systems can often be viewed as two coupled subsystems, one with oscillator-like degrees of 
freedom characterized by discrete frequencies and the other a wave-field characterized by continuous 
spectrum. When (artificially) decoupled from the wave-field, the discrete system has infinitely long- 
lived time-periodic bound states. Coupling leads to energy transfer from the system with oscillator 
like degrees of freedom to the wave field. In many situations, a (typically approximate) reduced 
description, which is a closed equation for the oscillator amplitudes, can be derived. This reduction 
captures the view of the oscillator degrees of freedom as an open system with an effective (radiation) 
damping term. In the problems considered in this paper, the reduced equation is of the simple form: 

id t A £ (t) ~ e 2 (A - iT )A e (t). (1.1) 

Here, e is a real-valued small parameter, measuring the degree of coupling between oscillator and 
field degrees of freedom. A e (t) denotes the slowly- varying complex envelope amplitude of the 
perturbed bound state. A is a real frequency and T > is an effective damping, governing the rate 
of transfer of energy from the oscillator to field degrees of freedom. 

For example, consider the general linear or nonlinear Schrodinger equation 

idt^ = H v <jf + e W(t,x,\(/) e \) <f, H v = -A + V(x). (1.2) 

Here, V(x) is a real- valued time-independent potential and W(t, x,\(p\) isa time-dependent potential 
(parametric forcing), W = (3(t,x), or nonlinear potential, e.g. W = ±|</>| 2 . Equation (1.2) defines 
an evolution which is unitary in L 2 (K). 

In this article we focus on the class of one-dimensional ionization problems, where 

W(t,x) = cos(fit) f3(x) 

where /3(x) is a spatially localized and real- valued function and /i > is a parametric forcing 
frequency. Thus, our equation is a parametrically forced Schrodinger equation: 

id t <J) e = H v <j) e + e cos(^) p(x) <jf. (1.3) 
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We focus on the case where the parameter e is real-valued and assumed sufficiently small. 

Assumptions for the unperturbed problem, e = 0: Initially we assume that the potential 
V{x), decays sufficiently rapidly as \x\ — > oo, although we shall later restrict to potentials with a 
fixed compact support. Furthermore, we assume that Hy has exactly one eigenvalue Ay < 0, with 
corresponding (bound state)eigenfunction, ipy(x): 

Hy lj)y = Ay Ipy, HVvlh = 1- (1-4) 

Thus, <f>°(x,t) = e- tXvt i/jy(x) is a time-periodic and spatially localized solution of the unperturbed 
linear Schrodinger equation: 

ldt<f> = Hy(f) 

We indicate an explicit dependence of Xy and tpy on V, since we shall be varying V. 

Fermi's Golden Rule: We cite consequences of the general theory of [33, 19, 20]. If //, the forcing 
frequency, is such that Ay + \i > 0, then for initial data, <j>(x, 0) = tpy(x) (or close to i/jy), the 
solution decays to zero as t — > oo. On a time scale of order e -2 the decay is controlled by (1.1), i.e. 

\A(t)\ ~ |A(0)| e- e2r ^\ 0<t< 0(e~ 2 ) (1.5) 

where \A(t)\ = \{i/jy, 4> e (t))\ and T[V] is a positive constant. Thus, we say the bound state has a 
lifetime of order (e 2 • T[V]) 1 and the perturbation ionizes the bound state. 

The emergent damping coefficient, T[V], is often called Fermi's Golden Rule [35], arising in the 
context of the spontaneous emission problem. However, the notion of effective radiation damping 
due to coupling of an oscillator to a field has a long history [22]. In general, T [V] is a sum of 
expressions of the form: 



(ey(-, Av es (Ay)),£ny(Vv))L2(Rd) = \ty(k res )\ 2 (fv(-,h-ea(^v)),Gw(i>v))L*(®.*) ) (I- 6 ) 

(see (3.2)) where Gwfyv) depends on the coupling perturbation W in (1.2) and the unperturbed 
bound state, ipy. Here, ey(-,k res ) = ty(k res )fy(-,k res ) is the distorted plane wave (continuum 
radiation mode) associated with the Schrodinger operator, Hy, at a resonant frequency k res = 
k res {\y), for which k 2 es G cr con t{Hy). ty(k) denotes the transmission coefficient and fy(x,k) a 
Jost solution. In Sees. 2 and 3 we present an outline of the background theory for scattering and 
the ionization problem, leading to (1.5), (1.6); see [33]. 

We study the problem of maximizing the lifetime of a metastable state, or equivalently, mini- 
mizing the scattering loss of a state due to radiation by appropriate deformation of the potential, 
V(x), within some admissible class, A\(a,b, fi): 

v min T[V\. (1.7) 

We refer to Eq. (1.7) as the potential design problem (PDP). Our admissible class, Ai(a,b, fj), is 
defined as follows: 

Definition 1.1. V G Ai(a, b, fi) if 
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1. V has support contained in the interval [—a, a], i.e. V = for \x\ > a 

2. V £ iJ x (K) and \\V\\ H i < b 

3. Hy = — d^ + V(x) has exactly one negative eigenvalue, Ay, with corresponding eigenfunction 
ipV £ L 2 (M), which satisfies Eq. (1.4) : Hyipy = Xyipy, [l^vlb = 1- 

4. k res = y/Xy + M > (formal coupling to continuous spectrum) 

Remark 1.1. Based on our numerical simulations, we conjecture that the hypothesis 2., imposing 
a bound of V, can be dropped. 

The idea of controlling the lifetime of states by varying the characteristics of a background 
potential goes back to the work of E. Purcell [29, 30], who reasoned that the lifetime of a state can 
be influenced by manipulating the set of states to which it can couple, and through which it can 
radiate. 

Remark 1.2. We discuss the potential design problem where 

1. (3(x) is a fixed function, chosen independently of V, for example, fi{x) = l[_2,2]( a; ) 

2. p[x) = V{x). 

Remark 1.3. How does one minimize an expression of the form (1.6)? 
We can think of two ways in which (1.6) can be made small: 

Mechanism (A) Find a potential in A\ for which the first factor in (1.6), \ty(k res )\ 2 is small, 
corresponding to low density of states near k 2 es . 

N.B. As proved in Proposition 2.6, |ty(A;)| > 0(e~ Ka ) for V with support contained in [—a, a]. 

Mechanism (B) Find a potential in Ai which may have significant density of states near k^ es (say 
\ty(k res )\ > 1/2) but such that the oscillations of fv(x, k res ) are tuned to make the matrix element 
expression (inner product) in (1.6) small due to cancellation in the integral. 

Indeed, we find that both mechanisms occur in our optimization study. 

Remark 1.4. We are interested in the problem of deforming V within an admissible set in such a 
way as to maximize the lifetime of decaying (metastable) state. Intuitively, there are two physical 
mechanisms with which one can confine wave-energy in a region: via the depth of the potential 
(material contrast) and via interference effects. We shall see that our (locally) optimal solutions, 
of types (A) and (B) find the proper balance of these mechanisms. 

1.1 Overview of results: 

1. In Proposition 4.11 we show that the optimal solution to Eq. (1.7) exists, for an admissible 
set, A±(a,b, n), derived from A\(a,b,fi) by relaxing a discrete constraint; see Eqs. (4.2) and 
(5.1). 
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Figure 1: Numerical demonstration for Eq. (1.3) (e = 1) of bound state time-decay for a typical 
potential well (left) and bound state persistence for an optimized potential (center). We plot 
the potentials Vi n u and V op t (black), corresponding ground states ipimt and ip op t (red), and forcing 
function j3{x) = lr_ 2 2]{ x ) (blue). The rightmost figure displays the time evolution of the projection 
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\{4>(t, •), tpv ("))| 2 f° r each potential. 
T[Vopt\ = 3.3 x 1(T 9 . 



Details are given in Sec. 6.6. V[Vi 
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2. Fix the admissible set A\(a, 6, fi), i.e. parameters a, b, fi, 5. Select an initial potential, Vinit £ 
Af(a, b, fi). For example, we have chosen a potential of the form V- in it = — Asech(Bx) li x i< 
with the parameters A and B appropriately chosen. We use a quasi-Newton method within 
A{(a,b,fi) and, after about 50-100 iterations, find a potential V op t G Af(a,b, fi), for which 
r[V] achieves a local minimum in Af(a, b, fi). 

3. In a typical search T[V ini t] ~ 10~ 2 and r[V opt ] ~ 10" 9 . Therefore, by Theorem 3.1 [33, 19, 20], 
the decay time for the solution of (1.3) with potential V = V op t and data (p e (0) = ipv op t is 
much, much longer than that for the Schrodinger equation with potential V = Vinit and 
data 4> £ (0) = ipv ini f Thus our optimization procedure finds a potential which supports a 
metastable state which has a very long lifetime, in the presence of parametric forcing coupling 
to scattering states. 

4. As an independent check on the performance of our optimal structures, we numerically solve 
the initial value problem for the time-dependent Schrodinger equation (1.3) with e = 1 for 
both V = Vina with data <f) e (0) = ipv init an d V = V opt with data ^ e (0) = Vv opi - Figure 1 
displays a representative comparison of these numerical experiments, revealing the decay of 
the bound state for Vi n u and a striking persistence (non-decay) of the bound state for V op t- 
The details of this simulation are given in Sec. 6.6. 

5. In Section 6.1, we investigate the optimization for classes of potentials with increasing sup- 
port, i.e. A{(a,b,fi) for an increasing sequence of a — values: < a\ < 0,2 < ■ ■ ■ < 
a m . Figure 2 shows local optima found in Af(aj,b, fi). As a is taken larger, the sequence 
V pt,a! , V op t A 2 ) • • • V op t t a m appears to take on the character of a truncation to the interval [—a, a] 
of a periodic structure with a localized defect. This suggests the following 

Conjecture 1.1. {V op t A } converges to V op t t oo(x) = V per (x) + Vi oc (x), where V per is periodic 
on R and Vi oc (x) is spatially localized. 

6. Our computational methods find locally optimal solutions which have small values of T due to 
either of the mechanisms discussed in Remark 1.3 above. In Sec. 6.3, we find the confinement 



5 



properties of potentials, which are optimal due to the cancellation mechanism ( mechanism 
(B) ), are very sensitive to perturbations in the forcing frequency away from the forcing 
frequency, fi, for which the optimization is carried out. 

7. In section 6.6 we study the stability or robustness of the state, ipv op ti f° r a locally optimal po- 
tential, Vopt- Time-dependent simulations of the parametrically forced Schrodinger equation 
are performed for an un-optimized potential, Vmit, and initial data ipv init + noise and for 
Vgpt, and initial data ipv opt + noise. Optimal structures effectively filter noise from a ground 
state, while a generic potential does not. The time scale of bound state radiation damping 
~ (e 2 r[Vopt]) is 3> the time scale for dispersion of noise. 

8. Our computations show that the inequality constraints of the (relaxed) admissible set, (4.2), 
are not active at optimal potentials. This is in contrast to studies of other spectral optimiza- 
tion problems, i.e. scattering resonances [13, 14] and band gaps [18, 32] and other problems 
[21, 7, 26, 11], where periodic patterns attaining material bounds are obtained. 

1.2 Outline of the article 

In section 2 we introduce the needed scattering theory background to explain resonant radiative 
time-decay and Fermi's golden rule, T[V], which characterizes the lifetime of metastable states. In 
section 3 theory we summarize the theory of [33, 19, 20] in the context of the ionization problem 
(1.3). In section 4 we introduce an appropriate regularization, Af, of the admissible set of po- 
tentials, A\ (see Definition 1.1) which is advantageous for numerical computation, and prove the 
existence of a minimizer within this class. In section 5 we outline the numerical methods used to 
solve the optimization problem. In section 6 we present numerical results for optimal structures 
and, as an independent check, investigate the effectiveness and robustness of these structures for the 
time-dependent Schrodinger equation with optimized and un-optimized potentials. Section 7 con- 
tains discussion and conclusions and Appendix A contains the detailed computations of functional 
derivatives and gradients used in the optimization. 

1.3 Notation and conventions 

1. L 2 (IR) inner product: (f,g) = J R f(x)g(x) dx 

2. l? co mp (M) is the space of compactly supported L 2 (M) functions and L 2 oc (IR) is the space of 
functions which are square-integrable on any compact subset of R. 

3. Weighted I? space: 

L 2 ' S (R) = {/: (1 + \x\ 2 )^f £ L 2 (R)}, s£R 
with norm ||/|| 2 2 , a = / R (l + |x| 2 ) s / 2 dx 

4. Weighted Sobolev space: 

H k ' s (R) = {/: d£f G L 2 ' S (M), < a < k}, seR 
with norm \\f\\% ktS = [| (1 + \x\ 2 )if\\ Hk 

5. B(X, Y) denotes the space of bounded linear operators from X to Y and B(X) = B(X, X). 
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6. Summation notation: ^_|_ f± = / + + /_ . 

7. The letter C shall denote a generic constant. 

Acknowledgements. B. Osting was supported in part by US NSF Grant No. DMS06-02235, 
EMSW21- RTG: Numerical Mathematics for Scientific Computing. M.I. Weinstein was supported 
in part by U.S. NSF Grants DMS-07-07850 and DMS-10-08855. The authors wish to thank P. 
Deift, T. Heinz, A. Millis, Y. Silberberg and C.W. Wong for stimulating conversations. 

2 Spectral theory for the one-dimensional Schrodinger operator 
with compact potential 

In this section, we discuss relevant properties of the Schrodinger operator Hy = —d 2 + V for 
sufficiently regular and compactly supported potentials, e.g. V G Ai(a, b, //). More general and 
complete treatments can be found, for example, in [1, 34]. 

2.1 The outgoing resolvent operator 

Let 7^ k G C. For V = 0, we introduce the outgoing free resolvent 

iP(x) = R [f](x,k)= [ G (x,y,k)f(y)dy, G (x, y, k) = i(2k)~ l exp(ik\x - y\) (2.1) 

JR 

defined for / G L^ omp (M). The function ip = Ro(k)f satisfies the free Schrodinger equation and 
outgoing boundary condition 

(-d 2 x - k 2 )iP = /, lim (d xT ik)ip = 0. 

x— >±oo 

For V / we introduce the outgoing resolvent, Ry(k) = (Hy — fe 2 ) -1 , satisfying 

(H v - k 2 ) oR v (k) = Id (2.2) 

and which, for Qk > 0, is bounded on L 2 (M) except for a discrete set of the form, ki = iK[, ki > 0, 
where — k? are eigenvalues of Hy. We have the identity 

R = R o (Hy - k 2 ) oR v = (Id + R V) o Ry, 

or equivalently, the Lipmann-Schwinger equation, 

R v = (ld + RoV)- 1 o Rq, %k>0, k^im (2.3) 

Proposition 2.1. The following are properties of the resolvent, Ry [1, 34]- 

1. The family of operators Ry(k): L 2 C0 (K) — > L 2 0C (M.), given by Eq. (2.3), exists and has a 
meromorphic extension to k G C. It has no pole for \ {0}. 

2. For k G M and arbitrary f G L 2 omp , the function ip = Ry(k)f satisfies 

(H v -k 2 ^ = f, lim (9,T^ = 0. 

x— >±oo 
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We denote by Gy(x, y, k), the Green's function, defined as the kernel of the integral operator Ry(k), 
in analogy with Eq. (2.1). In the upper half plane, Qk > 0, Gy(x,y,k) has a finite number of 
simple poles at ki = £kj, k; > . In the lower half plane, Qk < 0, Gy(x, y, k) may have poles at 
resonances, values of k for which the scattering resonance spectral problem: 

(H v - k 2 )i> = 0, lim (d x =F ik)ip = 0. 

x— >±oo 

has a non-trivial solution. 

A consequence of Theorem 4.2 of [1] is the following: 

Theorem 2.2 (Limiting absorption principle). For ^sk > 0, the resolvent Ry(k) = {Hy — k 2 )^ 1 
is a meromorphic function with values in B(L 2 ). For s > \, it can be extended to ^sk > as an 
operator on B(L 2,S , H 2 ~ s ) with limit 

Ry(k ) = lim Ry{k). (2.4) 

Sfc>0, 

fc->fc eK 

Throughout this paper, we shall understand Ry(ko), for ko S R, to be the limit taken in this way. 

Since we are interested in how the properties of solutions change with the potential, we make 
use of the resolvent identity 

Ry - R v = R V (U -V)R V . (2.5) 

We now refine Thm. 2.2 by showing that Ry(k) : L 2,s — > H 2 ~ s is (locally) Lipschitz continuous 
with respect to V. To prove this, we shall use the following bounds, used in the proof of Theorem 
2.2 [1]: 

\\Ro(k)\\ L 2, s ^ H 2,- s < C (2.6a) 
||(Id + i? (A;)y)- 1 || H 2,-^ H 2,- a < C(V), %k>0, s>* (2.6b) 

and the following 

Lemma 2.3. Suppose f € L 2 ~ S (R) has compact support with supp(f) C [—a, a]. Then f G L 2,S (IR) 
and 

\\f\\ L2 ,s<C(a)a 2s \\f\\ L2 ,-s. 
Proof. \\f\\ 2 L2 , s = J7 2 (l + \x\ 2 ) s dx < C(a)a 4s ff 2 (l + \x\ 2 y s dx = C(a)a 4s \\f\\ 2 L2 ,. s . □ 
Proposition 2.4. Fix a, b, \i G M, V € Ai(a, b, fj,) and for p > denote by 

B 00 (V lP ) = {U G Ai(o,6,m): \\V - U]^ < p}. (2.7) 
There exists a po > such that ifU€ B°°(V, po), then for s > \, 

\\R v (k) - Ru(k)\\ L 2,s^ H 2,s <C(V, Po ,a)\\V -UWco (2.8) 
uniformly for all k G R. 
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Proof. Let / G L 2 ' S (M). Using Eq. (2.5) and Thm. 2.2, we compute 

\\(R V - Ru)f\\ H2 ,-s < C(V)\\(U - V)Ruf\\ L 2,s. 

Then using Lemma 2.3 we have 

\\{U - V)Rjjf\\ L 2,s < C{a)a 2s \\U - V||oo|| J^^lli^.— 
<C{a)a 2s \\U -VW^WRufW^.s 

so that 

\\(R V - Rv)f\\ H 2,- s <C(V 1 a)\\U-V\\ x \\R u \\^H',-4f\\v».'' 
We now claim that there exists a po > and constants C(V) and C(V,a) such that for U G 

WRuWv^h*,- < C(V) / ■ (2.9) 

1 - C(V, a)po 

Equation (2.8) now follows once we have shown Eq. (2.9). To show Eq. (2.9), we use the resolvent 
identity 

Ru = (Id + (Id + RoV^RoiU - V)Y l R v 

and Thm. 2.2 to obtain 

WRuWv^h*.- < C(V)\\ (Id + (Id + RqV)- x Bq{U - V))' 1 11^,-^2,-. (2.10) 

Using Eqs. (2.6a) and (2.6b) we have 

\\(ld + RoV^RoiU -V)\\ H 2,s^ H 2,s < C(V,a) Po 

and Eq. (2.9) follows from using the Neumann series in Eq. (2.10). □ 

Proposition 2.5. Let V G A\(a,b, fi), k G M, k 0, s > \. There exists a po > such that if 
k'eB(k, Po ) 

\\Ro(k)- /2o(fc0llL2.»-*fla,-. <C(p ,a)\k-k'\ (2.11a) 

\\Rv(k) - Rv(k')\\ L 2, a ^ H 2,- s < C(p ,V,a)\k - k'\ (2.11b) 

Proof. Eq. (2.11a) follows from Eq. (2.1). To show Eq. (2.11b), we compute 

\\R v {k) - R v {k')\\ L 2, s ^ H 2,-s <|| [(Id + Ro(k)V)- 1 - (Id + Ro(k')V)- 1 ] R (k)\\ L 2,s^ H 2,s 

+ || (Id + R (k')V)~ 1 [Ro(k) - R (k')]\\ L 2,s^ H 2,s 
^CWild + Ro^V)- 1 - (Id + R (k')V)- l \\ H 2,~s^ H 2,s 
+ C(V)\\[Ro(k)-R (k')]\\ L 2,s_ H 2,-s (2.12) 

where we used Eq. (2.6). We now use the resolvent identity 

(Id + Ro(k)V)- 1 - (Id + Ro(k')V)- 1 = (Id + Roi^Vy 1 [R (k) - R (k')] V(ld + Ro(k')V)- 1 
and Eq. (2.6) on the first term in Eq. (2.12) to obtain 

||(Id + J R (A : )U)- 1 - {ld + Roik^Vy^s^.s < C(V)\\ (R (k) - R (k')) 
Now applying Eq. (2.11a) to Eq. (2.12) yields Eq. (2.11b) as desired. □ 
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2.2 Distorted plane waves, ey±(x]k), and Jost solutions, fy±(x; k) 

Distorted plane waves, ey±(x; k), are states which explicitly encode the scattering experiment of a 
plane wave incident on a potential resulting in reflected and transmitted waves. The Jost solutions, 
fv±( x i k)-, can be thought of as the states to which e ±lkx deform for nonzero V(x) in the spectral 
decomposition of Hy. In this section, we introduce these states and give their basic properties. 

The continuous spectrum of Hy is a c (Hy) = [0, oo). Corresponding to each point k 2 £ a c (Hy) 
are two distorted plane waves ey±(x,k) satisfying 

Hyey±(x,k) = k 2 ey±(x, k) (2.13a) 
lim (d x =F tk)ey(x, k) = 0. (2.13b) 

x— >±oo 

For V = these are the plane wave solutions eo±(x, k) = e ±tkx . For 7/0, the unique solution to 
Eq. (2.13) is given by 

ey±(x,k) = e ±lkx -Ry[Ve ±(;k)](x,k). (2.14) 

If V is compactly supported within [—a, a], for x ^ [—a, a], the solutions ey±(x,k) are given in 
terms of the transmission ty{k) and reflection ry{k) coefficients 



e %kx + ry(k)e lkx , x < —a 
ty(k)e tkx , x > a 



For k ^ 0, we have |ry(A;)| 2 + |iy(A;)| 2 = 1. If V is a symmetric, then ey-(x,k) = ey + (—x,k). 

The following proposition establishes that if V is compactly supported then |ty(fc)| is bounded 
away from zero, uniformly in k. We shall use this result in the interpretation of our numerical 
computations in Section 6.3. 

Proposition 2.6. Suppose supply) C [—a, a], k ^ 

IM*0I > exp f-min{l/|Jfe|, 2a} j |V(s)|daV (2.16) 
Proof. Consider the integral equation governing ey+(x, k): 

e v+ (x, k) = t v {ky kx - ^ Smk( - X k ~ y) V(y)e v+ (y, k) dy, x < a. 
For x > a, ey+(x, k) = ty(k)e lkx . Since k~ l sin(k(x — y)) is bounded by min{|A;| _1 , \x — y\} we have 
\e v+ (x, k)\ < \ty(k)\ + f mm{\k\-\ \x - y\} \V(y)\ \e v+ (y, k)\ dy (2.17) 
Therefore, by Gronwall's inequality 

\e v+ (x,k)\ < \t v (k)\ exp(f mm{\k\-\\x-y\} \V(y)\dy\ 

< \t v (k)\ exp(mm{\k\- l ,2a} \V{y)\d V y x < a, (2.18) 



10 



and thus 



> \e v+ (x,k)\exp (- minjlA;^ 1 , 2a} j \V(y)\ dy\ , x < a. (2.19) 



To bound \ev+(x, k)\, observe that for fixed k ^ 0, we can choose x* < —a such that arg(ry(fc)) = 
2kx* . Therefore 

\e v+ (x*,k)\ = \e lkx + r v {k)e~ lkx \ = |1 + r v {k)e~ 2lkx * \ = |l + |rv(fc)|| > 1- (2.20) 
The bounds (2.19) and (2.20) imply (2.16). □ 

The following proposition states that we can choose a constant to bound the distorted plane 
waves for all potentials in a small L°°-neighborhood of a V G .4.1 . 

Proposition 2.7. Fix a,b,p eR and V G A\(a,b,p) and let B°°(V,p) be as in Eq. (2.7). There 
exists a po > such that for U G B°°(V, po) the distorted plane waves eu±(x, k) satisfy 

l|e£/±(-,&)|U°°([-a,a]) < C{a,V,p ). 
Proof. Using Eq. (2.14), we compute 

„ikx-\\\ ^ /~</„M|/i i |„,|2\ — s r> \jTJikx 



\\Ru[Ue™ x ]\\ L ^ a>a]) < C(a)\\(l + \x\T S Ru[Ue^)\\ L ~ 
<C(a)\\(l + \x\ 2 r s Ru[Ue lkx )\\ H 2 
= C(a)\\R u [Ue ikx ]\\ H 2,-s 
<C(a,V, Po ) 



This last line follows from a Proposition 2.4. □ 

Definition 2.1. The Jost solutions, fv±(x,k), associated with the time-independent Schrodinger 
equation {Hy — k 2 )u = are defined by 

e v+ (x;k) = t v (k) f v+ (x;k), e V -{x;k) = t v (k) f v _(x;k), (2.21) 

where fy+(x; k) ~ e lkx as x — > +oo and fv-(x; k) ~ e ~ lkx as x _ OCm 

By results of [8], for any k G M and any compact subset, C, of M 

max|/y±(s;fc)| < K KC < oo (2.22) 

Note also that Propositions 2.7 and 2.6 imply a bound on \fv±\ in the case where V has compact 
support. 
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2.3 Spectral decomposition of the ID Schrodinger operator 

We state the spectral theorem in terms of the distorted plane waves (see e.g. [34]): 

Proposition 2.8 (Spectral Decomposition). Let ey± and fy± denote the distorted plane waves 
and Jost solutions given by (2.14) and (2.21). Let Xj for j = 1.. .N be the eigenvalues of Hy with 
corresponding (normalized) eigenf unctions ipj(x). Then, h = P^h + P c h where Pd and P c are, 
respectively, projections onto the discrete and continuous spectral parts of Hy given by 

1 f°° 

P c h=— { (e v+ (-,k),h)e v+ (x,k) + (e V -(-,k),h)e V -(x,k) } dk 
1 f°° 

= 7T l(fv + (;k),h)fv + (x,k) + (fv„(;k),h)fv„(x,k)] \t v (k)\ 2 dk (2.23) 
^ Jo 

N 

Pdh = ^2>^j(tpj, h)ilfj(x) 

3=1 



Moreover, 



1 



oo 



g(Hy)h=— J 9 {k 2 ) [ (f V+ (;k),h) f V+ (x,k) + (f V -(;k),h) f V -(x,k) ] \t V (k)\ 2 dk 

N 

+ X>( A i)(Vy.fr)^0»0 > ( 2 - 24 ) 

3=1 

where g is any B or el function. Finally, by approximation we have that (2.24) holds with = 5((), 
the Dirac delta distribution in the distributional sense. 



3 Radiation damping and Fermi's Golden Rule 

In this section, we explain how r[V], given in Eq. (1.6), emerges as the key quantity controlling 
the lifetime of the metastable state. We now state a theorem on the ionization and decay of the 
bound state and then sketch the idea of a proof, which explains the mechanism of decay and (1.5). 
A detailed proof can be found in [33, 19, 20]. The following result holds for generic potentials with 
one bound state. In particular, these hypotheses are satisfied by V G Af(a, b, fx). 

Theorem 3.1. Consider the parametrically forced Schrodinger equation 

id t( j) e = H V (p e + e cos(fj,t) p[x) (j) e . (3.1) 

Assume V and (3 satisfies the general conditions of [33, 19]. Consider the initial value problem for 
Eq. (1.3) with <j) € (x,0) = (f> € L 2 ' a (R), where a > 1. A ssume 

1. k 2 ? = Ay + fi > (resonance with the continuum at 0(e 2 )) 
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T[V] = j (Mv, 5(H V - kl)P c Pi>) (3.2a) 
= 16 k v 



2. T[V] > 0, where T\V] is the non-negative quantity defined by 

(^y,5(H v -k 2 v )P c ^) 

l ] -Y,\(^V,ev±(;k v ))\ 2 (3.2b) 

= WW |ty(%)|2 El^'^±(-,A; y ))| 2 , (3.2c) 

where ey± and fy± denote, respectively, the distorted plane wave and Jost solutions, and 
ty(k) denotes the transmission coefficient. 

Then, there exists eo > such that for e < eo 

\{i, v A e {:t))\-\^v,^)\e-^ [V]t + 0(e), 0<t<O(e~ 2 ) 

\\<F(;t)\\v,- < rl||^|| i2 , CT , t»i. 

Remark 3.1. For certain choices of potentials, the parametrically forced Schrodinger (ionization) 
problem is exactly solvable by Laplace transform methods and the time-behavior can be computed 
for all e > 0. See, for example, [6, 5]. 

A sketch of the proof. In this sketch, we drop the subscript on ipy and superscript on <j) e . For 
small e, it is natural to decompose the solution as 

4>{t,x) = a(t)ip(x) + cf> c (t,x) (3.3) 

where a(t) = (tp, (/)(■, t)) and cp c = P c [4>] is the continuum projection; see (2.23). To simplify the 
discussion we take as initial data: 

a(0) = a , M0,x) = 0. (3.4) 

Substitution of (3.3) into (3.1) and projecting onto the discrete and continuous spectral parts of 
Hy yields the following coupled system: 

(idt — X)a(t) = e cos(/ii)(V>, (3ip)a(t) + e cos(/j,t) (tp , [34> c ) (3.5a) 
(tdt - H v )<t) c = ecos(fit)P c [^}a(t) +€COs(/it)P c ^ c ]- (3.5b) 

Since e has been assumed small, the coupling between a(t) and 4> c (t,x) is weak. We now proceed 
to make a set of simplifications leading to a minimal model, in which the mechanism of radiation 
damping is fairly transparent. First, since the first term on the right hand side of (3.5a) contributes 
an order e mean-zero frequency shift from A, we neglect it. Second, from equation (3.5b) we formally 
have that <p c = 0(e). Therefore, the last term on the right hand side of (3.5b) is C(e 2 ) and we 
therefore neglect it. Finally, the second equation evolves in the continuous spectral part of Hy and 
we formally replace Hy by Hq = — A. 

The resulting system is the following Hamiltonian system of an oscillator of complex amplitude 
a(t) coupled to a field 4> c (t, x): 

{idt — X)a(t) = e cosd^t) (ip , /3(f> c ) (3.6a) 
(tdt + A)0 C = e cos(nt)pipa(t). (3.6b) 
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We can exploit a separation of time-scales by extracting the fast phase from a(t) via the substitution 

a(t) = e~ tXt A(t), 

giving the following equation for the slowly varying amplitude, A(t): 

id t A{t) = ecos(fit)(i;,f3<j) c )e lXt (3.7) 
Now, Duhamel's formula is used to rewrite Eq. (3.6b) as 

<j) c (t) = -ie I e jA( *- s) cos(/is)/3V>a(s)ds 
Jo 

since </> c (0) = 0. We insert this back into Eq. (3.7) to obtain the closed equation for A(t). 

8 t A{t) = -e 2 cos(fit)e- lXt A(t) [ (f3ip,e lA ^- s ^tp)cos(fis)e- lXs A(s)ds 



Writing cos([it) = \ (e lfJ-t + e JAt *), we find that if k 2 es = A + /i > 0, then it is a resonant frequency 
and 

d t A(t)^--e 2 e- tk ^ t A(t) / (J3i/}, e lA{t ~ s) P^e' 1 ^ 3 A{s) ds 

« --^(^^-A-kl^-ity-^^Ait) 

Here, (—A — E — iO)" 1 = lim,5|o( — A — E 2 — id) -1 . The choice of regularization is dictated by the 
outgoing radiation condition for t — > +oo; see [33, 19]. 

Returning to the original (un-approximated) equations (3.5), we have analogously 

8 t A(t) « -\e 2 {^, {H v - k 2 res - i0y l P c [^])A{t) = -e 2 (A + %T)A(t). (3.8) 

The coefficient of A(t) in (3.8) can be computed by applying the functional calculus identity (2.24) 
to the function g(s) = (s — k 2 . es — it) -1 , together with the distributional identity 

lim(s - k 2 es - it) -1 = P.V. (s - k 2 ^)' 1 + ^7^ 5(s - k 2 es ) 
and the identification s — > Hy. In particular, 

= \-^( Mv,6(H v - k 2 res ) P C M > 

1 f°° r i 

= ^ / t(k 2 -k 2 res ) [\(f v+ (.,k),^ v )\ 2 + \(fv-(;k),Wv)\ 2 \ \t v (k)\ 2 dk, 

from which the expression (3.2) follows after setting v = k 2 and carrying out the integral. 
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4 A constrained optimization problem: design of a potential to 
minimize radiative loss 

We now consider the Potential Design Problem (PDP) given in Eq. (1.7) with T[V] defined in Eq. 
(3.2). We begin by discussing the set of admissible potentials A±(a, b, fi) defined in Def. 1.1. For the 
purpose of numerical computation we relax he admissible set, A\ — > A\ , by replacing the discrete 
constraint (Hy has exactly one eigenvalue) by an inequality constraint, in terms of a regularization 
parameter, 5. We then show that the objective function is locally Lipschitz and that a solution to 
the PDP exists in the modified admissible set, Af(a,b, /i). 



4.1 The admissible set A\ and its relaxation, A{ 

Denote the Wronskian of the distorted plane waves, ey±, by 

W v (k) =Wron(e v+ {;k),e V -(;k)), keC. (4.1) 

Zeros of Wy{k) correspond to poles of the Green's function Gy (x, y, k) as introduced in Sec. 2.1. In 
particular, the zeros of Wy(k), in the upper half plane, are eigenvalues. The number of eigenvalues 
is increased or decreased by one, typically through the crossing of a simple zero of Wy{k) through 
k = as V varies. 1 Our strategy to fix the number of eigenvalues is then to start with a one bound 
state potential and deform V, keeping Wy(0) 7^ 0. However, numerically it is advantageous to 
replace contraint Wy(0) / by the inequality constraint Wy(0) 2 > 5. For 5 > 0, we regularize 
A\(a,b, fi) by introducing 

A^a, b, n) = Ai(a, b, fx) D {V: Wy{0) 2 > 5}. (4.2) 

Remark 4.1. Note that the set of admissible potentials A\ (a, b, fi) is not convex. Indeed, counter- 
examples can be explicitly generated and are illustrated by the following cartoon superposition of 
potential wells at sufficiently separated points [12]: 



0.5 



+ 0.5 



The two potentials on the left hand side of the equation support a single bound state while the 
convex combination on the right supports two. 



Lemma 4.2. For 5 > 0, let 7: [0,1] — >• Af(a,b, //) be a smooth function valued path. If 
supports a single bound state then so does H^^) f or t G [0, 1] • 

Proof. We show that no bound states are lost along the path j(t). A similar argument shows that 
no bound states are gained. 

Define f(k,t) = W 7 ( t )(fc) and consider the equation f(k,t) = 0. Let Ao denote the eigenvalue 

. 7^ 0, i.e. An is a simple eigenvalue, by the 



of iI 7 (o), i.e. f(iy/\Xo\,0) = 0. Since d k f(k,t] 



1 Potentials, V, for which Wv(k = 0) = are called exceptional. The value k — 0, corresponding to edge of the 
continuous spectrum is then called a zero energy resonance or a half-eigenvalue with half-bound state ev±(x, 0) [31]. 
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implicit function theorem, there exists T > such that for \t\ < T, there is a parameterized family 
1 1-)- At with At < and f(iy/\\t\, t) = 0. Let 



t* = sup{0 < t < 1: /(«\/|At|,t) = and A t < 0} > 0. 
If A t # = 0, t# < 1 then f{0,t*) = W (t#) (0) = 0. This contradicts 7(t) C Xf(a,6,/i) Therefore, 



i# = l. 

Let rjy(x) = ey±(aj, 0) be the distorted plane waves at A; = which satisfy 

#yf?y = 0, lim 7?y = 1; 



□ 



(4.3) 



see equation (2.15). 

Our gradient-based optimization approach requires that we compute the variation of the Wron- 
skian, WV(0), with respect to the potential V. This calculation will also be used to establish 
Lipschitz continuity of WV(0). 

Proposition 4.3. Let Tj v (x) satisfy Eq. (4-3). The Frechet derivative of the Wronskian Wv(0) = 



Wron(r]y , r] v ) : L 



camp 



with respect to the potential is given by 
5W V (0) 



5V 



-VvVv 



Proof. See Appendix A. 



□ 



Remark 4.4. If V is symmetric, then SW ^y°^ is symmetric. 

To prove that WV(0) is locally Lipschitz, we use the following lemma. 

Lemma 4.5. Let f[V] : L 2 ([— a, a]) — > M be a Frechet differentiable functional with 



U -V) +o(\\U -V\ 



2)- 



Suppose further that the variation |^ is bounded in an L°° -neighborhood ofV. Then there exists a 
po > and a constant C(po, V, a) such that for U G B°°(V, po) 

\f[U] - f[V}\ < C(po,V,a)\\U - V\\ Lao{ha>a]) 

Proof. The Mean Value Theorem and Proposition 4.3 imply that there exists a po > such that 
for every U £ B°°(V, po), there exists a potential V = tV + (1 — t)U for some t £ [0, 1] such that 



,u-v 



This gives the estimate 



r o~f 


dx 




V 



\f[U)-f[V]\<\\U-V 
The proof is completed by choosing C(po, V, a) = sup V( - Boo ^ v po ^ 



a sf 



r 

J— a 



5V\V 



dx 



□ 
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Proposition 4.6 (local Lipschitz continuity of Wy(0)). Fix a,b,p G R, 5 > 0, and let V G 

Af(a,b, p). For p > denote by 



B°°(V,p) = {U £ A{(a,b, p): \\V - < p) 



(4.4) 



the L°°(1R) ball around V in Af(a,b,p). Let Wy(0) be as defined in Eq. (4-1)- There exists po > 
and a constant C(po,V,a) such that if U G B°°(V,po) then 

\Wu(0) - W v (0)\ < C( Po ,V,a)\\U- V\\oo. 



SWy 



Proof. Propositions 2.7 and 4.3 give that -m 
of V. The result now follows immediately from Lemma 4.5. 



r]yT]y is pointwise bounded in a neighborhood 

□ 



4.2 Properties of the objective functional, r[V] 

In this section, we begin with a formal calculation of the Frechet derivative of r[V], given by (3.2). 
We then show that r[V] is (locally) Lipschitz with respect to V. 



Proposition 4.7. Let L > a. The Frechet derivative of T\V] \ L 2 comp (\— a, a]) 
(3.2) with respect to the potential V is given by 



SV 



5T 

1 



[5V] 



Slpy 



^ v RvWX v )P c 



5T 

Se v ± 



5ey±(-,ky) 

SV 



[SV] 



8k v 



'R^2(ev±,(3^v}/3ev± 



r / Skv AT/\ 

kv { W> 5V) 



given in Eq. 

(4.5a) 
(4.5b) 



2k v 



((3i/j v , A±)*p v - ev±Rv{kv) Wv] 



where 



A±(x) = ± ixe lkvx - Ry(ky) 2k V (/) V ± ± ixVe ±lkvx 



+ 



lkyL 

2k v 



L >v±{-L)e v+ {x,k v ) + 4> v ±(L)e V -(x,k v )) 



(4.6) 



ey± = ey±(-, ky), and <ftv±(x) = e ±ikvX — ey±(x, ky) satisfies Eq. (A. 6). 

Proof. Equation (4.5a) is obtained by the chain rule. A detailed computation of each term is given 
in Appendix A. □ 

Remark 4.8. If the potential and j3 are symmetric, so is |R. 

Proposition 4.9 (local Lipschitz continuity of V ). Fix a,b, p. G M, S > and V G Af(a,b,p) 
and define B°°(V,p) as in Eq. (4-4)- There exists po > and a constant C(V, po,a) such that if 
U eB°°(V,p) then 

\T[U]-T[V}\ <C(V,p ,a)||E7--V|| oo . 
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Proof. First we use the triangle inequality to obtain 

|r[D] - F[V}\ = | -^\Wi>v,ev±(;ku))? - 1 ± z \(04'v,ev±(; k v ))\' 



~ 16 4^ 



1 1 

kjj ky 



16 V 

K/3^,e C /±(-,%))| 2 



+ 



+ 



^ E [ iK^.et/iCfct/^p-K^v.et/iO,^))! 1 



16A;y 



IWy,et/±(',fe)>| 2 - |(/3W,ev±(-,fcy))| 2 
+ |KWy,ey±(-,%))| 2 - |(^y,ev±(-,A;y))| 2 



(4.7) 



We now treat the terms A± — D± in Eq. (4.7) in turn. We'll repeatedly use the inequality 
IN 2 - |6| 2 | < \a + b\\a-b\. 



A. We compute 



\k v - ku\ 



lGkyku 

Recalling ky = \/^v + M an d using Eq. (A. 3) we have 

1 



\(Wu,eu±{;ku))\ 2 . 



\k v - k v \ <—\Xy- Xu\ + o(\\ v - Xu\) 

Iky 

< J-(^,|^-^|>+o(l|^-^[|oo) 



< 



2ky 
1 

2Av 



\U -VlU + oiWU -v\ 



B. We compute 



B± 



1 



< 



< 



16% 
1 

16Av 
1 

16% 



\{Mu,eu±(;ku))\ 2 - \(Mv,eu±(;ku))\'' 



\(P(ifo + ifr),e u± (-,lto))\\{Ptyu- , il>v),e u± (-,k u ))\ 
\{PO>U + 4>v), eu±(-, ku))\\\eu±(- , ku)\\L2(K)\\P(4>u - ipv)h- 



Now using Eq. (A. 4) and Thm. 2.2 we obtain 



'm-ipvh < \\R v (y/X)Pc[MU-V)]\\ 2 + o(\\U-V\ 
<C(V)\\U-V\\oo + o(\\U-V\\)oo. 



(4. 
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C. We compute 



C± 



1 



< 



1 



\W^v,e u± (;ku))\ 2 - mv,e v± (-,ku))\' 



\(0ipv,eu±(-,ku) + ev±(-,ku))\\WiJv,eu±(-,ku) - e v ±(-,kv))\ 



16 V 

Using Eq. (2.14) we have that 

eu±(;ku)-e v± (;ku) = -R u (k u )P^ k ° IB ] + Rv(ku)[Ve ±tk " x ] 

= Ru(ku)[(V - U)e ±lk ^\ + (Rv(ku) - Ru{ku))[Ve ±lk ^ x ] 

Fact: For K C -6(0, r) compact we have 



L°°{K) 



;i + |x| 2 )~2/(i + |x| 



\L° a (K) 



< C r \\(l + \x\ 



<C r \\(l + \x\ 2 )-*f\\m ( 
= C r \\f\\m -»(E)- 



Thus, by Prop. (2.4) and Eq. (4.9) we have 

\\eu±(-,ku) - ev±(-,ku)\\L°°(K) 

< C r \\e u ±(-,k u ) - e v ±(-,ku)\\ H i,- 

< C r (\\Rv\\l^h*.-4(V - U)e ±tk » x \\v.° + H^v-^lli^^-aliye^l 
<C r (\\R v \\ L ^ m ,-s\\e ±ik n\^HK)+C(V,^ 



D. We compute 



At 



l 



16% 
1 



|</3Vv, ev±(; k v ))\ 2 - \((3ip v , e v ±(; k v ))f 



< Y^;K/ 3 V ; V,ey±(-,%) + e v ±(;kv))\\(Pih , ,ev±(',ku) ~ e v ±{-,k v ))\ 



But, 



(4.9) 



e V ±(; ku) - e v ±(; k v ) = e ±lkvX ( e ±*(%-M* _ ^ _ [ Ve ±*k v x ^ e ±t{k v -k v )x _ 1 

+ (R v (k v )-R v (k u ))[Ve ± ^} 



(4.10) 



Insertion of (4.10) into the above bound on D± and use of Prop. 2.5, (4.8) and Proposition 2.7 
implies \D±\ < CWV-UWoo. 



Proposition 4.9 now follows from assembling the estimates A-D. 



□ 
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Remark 4.10. As mentioned in Remark 1.2, we shall consider optimization problems where (i) 
/3(x) is a fixed specified function and where (ii) /3(x) = V(x). For type (ii) problems r[V] is given 
by 

rm = I ^El(^y ) ey±(-,M}| 2 - (4-11) 

That r[V] in this case is Lipschitz follows by the same arguments as above. Furthermore, it is 
Frechet differentiable and the additional contribution of the Frechet derivative of r[V], to be added 
to the expression in Proposition 4.7, is given by: 



— i> v (x) 3? ^2 ev ± ( x ' k v) ( e v± {;kv), Vip v ) ■ 



8k v ± 
4.3 Existence of a minimizer 

We show that the potential design problem attains a minimum in the admissible class Af(a, b, //). 
Define 

7 f (a, b, n) = M{F[V] : V G A{ (a, b, fj,)} > 0. (4.12) 
Proposition 4.11. There exists V* € Af(a,b,fj.) such that T[V*] = 7^(0,6,//). 

Proof. Let {V n } C Af(a,b,/j,) be a minimizing sequence, i.e. lim^oo T[V n ] = 7*. Since ||V^||#i(R) < 
b, there is a weakly convergent subsequence converging to V* G H . Moreover, the family {V n } 
is uniformly bounded and equicontinuous on [—a, a]. By the Arzela-Ascoli theorem, there is a 
subsequence (which we continue to denote by {V^}) converging to V* G H 1 and such that V n — > V* 
uniformly on [—a, a]. By Prop. 4.6, Wy(0) is continuous with respect to V on A±(a, b, fj,), implying 
K, G Af(a, b, fx). By Prop. 4.9, r[V] is continuous on Af(a, b, fx) and 



r[K] = iimr[K] = 7* 



□ 



Remark 4.12. Since ^4* is not convext (Remark 4.1), uniqueness of the minimizer is not guaran- 
teed. 

Corollary 4.13. By Remark 1^.10, a minimizer also exists if we take T with (5 = V. 

5 Numerical solution of the optimization problem 

In this section, we discuss a numerical solution of the potential design problem (PDP) 

min T[V] (5.1) 

VeA({a,b,fi) 

for fixed a,b,/i,S, and /3(x), where T[V] is given in Eq. (3.2) and A((a,b,fj,) in Eq. (4.2). 
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Forward Problem. We refer to the evaluation of T\V] and Wy (0) for a given as 
the forward problem. 

To evaluate the objective function r[V], first the eigenpair (Ay, tpv) satisfying Eq. (1.4) is 
computed using a three-point finite difference discretization and the Matlab eigs command. Next, 
the distorted plane waves ey±(x, vV + Ay) are computed using the decomposition as in Eq. (A. 5) 
and then solving Eq. (A. 6) using the same discretization. The integrals for T[V] are then evaluated 
using the trapezoidal rule. 

The evaluation of Wy (0) requires the distorted plane waves at k = 0, which are computed using 
a Crank-Nicholson method. For numerical stability, the Wronskian, which is analytically a constant 
in x, is computed on a uniform grid and is averaged over the spatial domain. As a check on the 
discretization, we ensure that the variance of the Wronskian does not exceed a specified tolerance. 

Optimization Problem. Local optima of Eq. (5.1) are found using a line-search based L-BFGS 
quasi-Newton interior-point method [25] as implemented in the Matlab command fmincon. We use 
the optimize-then- discretize approach, where gradients are computed as in Proposition 4.3 and 4.7 
and evaluated using the discretized counterparts. The constraints, 



are enforced using a logarithmic-barrier function. The method terminates when the line search 
cannot find a sufficient decrease in the objective function. 

In the numerical experiments, presented in Section 6, we use a computational domain larger 
than the interval [—a, a] defining the support of the potential V and depending on the magnitude 
of a, between 1000 and 3000 grid points. The method converges in less than 100 iterations and 
takes approximately 5-20 minutes using a 2.4 GHz dual processor machine with 2GB memory. 

Time-dependent simulations. In Sections 6.6 and 6.7 we study time evolution for the initial 
value problem in Eq. (1.3). This is accomplished using the same discretization as above and the 
time stepping routine for stiff ordinary differential equations implemented in odel5s in Matlab. 
The outgoing boundary conditions are approximated using a large domain with a dissipative term 
localized at the boundary. 

6 Results of numerical experiments 

In this section, we present the results of many numerical experiments using the methods described in 
Section 5 to study locally optimal solutions of the potential design problem (5.1). The constraints 
in Eq. (5.1) depend on fi (forcing frequency), a (support width), b (H 1 bound on V), and 8 
(relaxation parameter), while the objective function depends on the choice of spatial perturbation 
of the potential, /3(x). For 6 sufficiently small and b sufficiently large, we find that in all numerically 
computed solutions of Eq. (5.1), a local optimum is achieved at an interior point of the constraint 
set, Af(a,b, fj,), i.e. the constraints (5.2) are not active at the optimal solution. This is in contrast 
to the structure of optimal solutions of other design problems studied in [7, 13, 14, 18, 21, 26] where 
the optimal potentials always attain the bounds and are referred to as "bang bang" controls. 
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Figure 2: Locally optimal potentials for varying values of support a, with fixed frequency [i = 2 
(left) and varying forcing frequency /i, with fixed support a = 80 (right). The potentials are 
symmetric in x, only x > is plotted. 



We conjecture that the H 1 bound on V can be relaxed in Proposition 4.11 and that the con- 
straints of a compactly supported potential with a finite number of bound states is sufficient for 
the minimization to be well posed, i.e. there exists bo, So > such that for b > bo and 5 < 5o- 

min r[V] = min r[V]. 
VeA{(a,b,n) VeAl(a,oo tt j,) 

Thus, we consider potential optimization problems for the two classes of /3(x) in Remark 1.2 and 
vary [i and a. 



6.1 Optimal potentials for varying support size, a, 
with forcing frequency /i = 2 and 0(x) = 1[— 2,2] (^) 

In Fig. 2 (left) we plot locally optimal potentials for 5 different values of the support, a. (The 
potentials for different values of a are shifted vertically.) The potentials that emerge are symmetric 
in x (see remarks 4.4 and 4.8) and periodic on the interval [—a, a] with a defect at the origin. 
Let V* denote the optimal potential for support parameter a. We note that existing structure 
changes very little as we increase a. That is, V* and V£ are nearly equivalent on the set [c, c] where 
c = min(a, b). This is numerical support for Conjecture 1.1. 

The following table gives the value of r[V] for the sequence of potentials in Fig. 2(left). 



a 


4 8 


16 


32 64 


r 


3 x 10~ iU 2 x 10" n 


1 x 10~ 9 


3 x 10" y 8 x 10~ iU 



The non-monotonicity of V with increasing support parameter, a, reflects the fact that we are only 
able to compute local minima of V. For small a, a relatively small number of design variables give 
good numerical accuracy in evaluating the oscillatory integral that defines T and the optimization 
method converges easily. However, as a increases, the numerical method degrades as we are forced 
to balance accuracy with the number of optimization variables. 
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6.2 Optimal potentials for varying forcing frequency, //, 
with fixed support size, a = 80, and = lf-2,2] 

In Fig. 2 (right) we plot locally optimal potentials for 5 different values of forcing frequency /i. The 
optimal potentials vary smoothly as we change /i with the period of the oscillation in the tails of 
the potentials decreasing with increasing /i. 

The following table gives the value of r[V] for the potentials in Fig. 2(right). 





1.6 


1.8 


2 


2.2 


2.4 


r 


2 X 10~ 9 


7 x KT 9 


2 x KT 8 


4 x i(r 8 


2 x io -8 



6.3 Two mechanisms for potentials attaining small T 

The functional to be minimized, r[V], is given by 

r[V} = ^j- \Mk v )\ 2 J2WvJv±(;k v ))\ 2 ; 

see (3.2). As discussed in the introduction, two possible mechanisms can be used to decrease the 
values of T[V]; see Remark 1.3: 

Mechanism (A) Find a potential in A\ for which the first factor in (1.6), \ty(k res )\ 2 is small, 
corresponding to low density of states near k 2 es , or 

Mechanism (B) Find a potential in A\ which may have significant density of states near k 2 es (say 
\ty(k res )\ > 1/2) but such that the oscillations of fv(x, k res ) are tuned to make the matrix element 
expression (inner product) in (1.6) small due to cancellation in the integral. 

In Fig. 3 we display the results of numerical simulations illustrating examples of both mecha- 
nisms at work. On the left is the potential, VA,opt(x), and diagnostics exhibiting mechanism (A) 
and on the right we exhibit mechanism (B) for the potential labeled VB )0 pt{x)- For both examples 
we choose j3(x) = l\-2,2]( x )- 

The potential VA,opt(x) is obtained via optimization on the set A\ with a = 64 and /j,a = 2 
(same as in Fig. 2 (left)). The potential VB : opt(x) is obtained via optimization on the set A\ with 
= 8 and \xb = 4. The first row of figures displays the potentials as functions of x. The value of 
r for V A , op t{x) and V B , op t(x) are T [V A ,opt] = 1.2 x 1(T 8 and T [V B>opt ] = 1.3 x 10" 15 . 

The second row of plots is of the transmission coefficients |t\/(fc)| 2 (see Eq. (2.15)) of VA,opt{x) 
and VB,opt( x )- The small vertical arrows along the k— axes indicate the location of the resonant 
frequency k = ky. 

Remark 6.1. Relevance of the transmission coefficient, ty(k), to the density of states: 

Consider a periodic potential, q(x), defined on R. The spectrum of — d 2 +q(x) is equal to the union of 
closed intervals (bands) separated by open intervals (gaps). Now consider qht{x) = q{x)l^M,M]{x). 
The decaying potential q\j{x) has continuous spectrum extending from zero to infinity. We expect 
however the spectral measure, associated with the self-adjoint operator, — d 2 + qM{x) for M > 1, to 
have little mass on those intervals corresponding to the gaps in the spectrum of the limit operator, 
—d 2 + q(x). Related to this is the observation that the t g i, UM Ak), for —d 2 + q(x)l[_M,M]( x ), is 
uniformly small, for k 2 in the spectral gaps of the limit operator, and converge weakly to one for k 2 
in the spectral bands; see, for example, [2, 17]. Thus, by plotting the amplitude of the transmission 
coefficient for our optimal potentials we can anticipate whether the density of states is small and a 



23 



Low DOS Mechanism- VA,opt Cancellation Mechanism- VB,opt 




Figure 3: Comparison of two locally optimal potentials achieving small T[V] due the low density 
of states mechanism (left) and the cancellation mechanism (right). First row: displays plots of 
potentials. Second row: transmission, ty(k), a measure of the density of states. At the resonant 
frequency, ky, indicated by the arrow, ty(k) is very small on the left and approximately one on the 
right. The left figure inset shows that the resonant frequency is distinct from the resonant spike 
in the "gap". Third row: distorted plane waves, \ey + (x, ky)\. Fourth row: r[V^,t;/x] vs. /x for 
V = V op t t A optimized for forcing frequency fiA = 2 (left), and for V = V op t,B optimized for forcing 
frequency /is = 4 (right). Note contrasting sensitivity to perturbations in fx away from (j-a,b- 
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spectral gap is being opened around the resonant frequency, ky. Thus, if ky lies in an interval of 
very low transmission, ty(k), the T, given by (3.2) will be small. 

The left plot in the second row shows that the transmission coefficient for VA,opt(x) is very 
close to zero very near the resonant frequency, ky Aopt . On the right we see that for VB,opt(x) the 
transmission coefficient very near ky A opt close to one. 

In the third row of plots, for each potential, we plot the modulus of the distorted plane wave 
at the resonant frequency, \ey + (x, ky)\. (Recall that for a symmetric potential, ey_(x,/c) = 
ey+(—x,k).) The modulus of the distorted plane wave associated with VA,opt{x) decays rapidly 
as it enters the support of the potential, as expected since the transmission coefficient is nearly 
zero for this frequency (see Eq. (2.15)). The modulus of the distorted plane wave associated with 
VB,opt(x) is nearly unity over the support of the potential. 

In the bottom row of plots of Fig. 3 we highlight an additional distinction between these two 
mechanisms. We fix the optimal potentials, v~A,opt{x) arid VB,opt(x), respectively optimized for 
forcing fixed frequencies ha and hb- We then consider the variation of the function fj, \— > T[V op t', 
where 



T[Vop t ;fj] = ^-j=== t Vopt (y/\ Vopt +n ) | ^ | ((3ip Vopt , f Vopt ± yj \ Vovt + 



Here, \x varies over a range of forcing frequencies above and below (j,a, respectively, fis- 

We find that for VA,opt(x), the value of T is relatively insensitive to small changes in /i near 
[i A- Indeed, this is expected. Small variations in //, imply small variations in y/Xy + fi. Therefore, 
if ky Aopt = sJVA,opt + fJ-A is located in a spectral "gap", then for values of \i near fj,A, k(fi) = 
V VA,opt + ft is also in this "gap" . Therefore, ty A opt (k(fJ,)) and therefore r[VA i0p t,//] is small. 

In contrast, for VB,opt{x), the range of [i for which r[Ve,opt; fA remains small is extremely narrow; 
the smallness of the oscillatory integral, T[VB,opt', m]j is not preserved over a range of values of fx. 

Remark 6.2. These observations on the sensitivity of T[V op t,fj] with respect to the forcing fre- 
quency, fi, for the two different kinds of optimizers, A— type and B— type, should have ramifications 
for applications. 

By Proposition 2.6, 

VeA{(a,b,n) => \ty(k)\ > exp (-4a 2 6) . (6.1) 

Thus we find that T > due to Mechanism (A). However, in principle, one could find a potential 
such that due to perfect cancellation, T = by mechanism (B). Indeed, the potential Vb has an 
extremely small value of T. 

6.4 Further discussion of mechanism (A); potentials which open a gap in the 
spectrum 

We have observed that some locally optimal potentials, e.g. the potential associated with the left 
column of Fig. 3, have small values of T due to mechanism (A), creating a low density of states at 
the resonant frequency ky opt . We explore this phenomena further here and discuss the relation to 
Bragg resonance. 
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Figure 4: (Left) For the sequence of potentials in Fig. 2 (left), a spectral gap forms as a f oo. 
(Right) For each of the potentials in Fig. 2(right), the resonant frequency lies in a spectral gap. 

For the sequence of potentials given in Fig. 2 (left) corresponding to an increasing sequence 
of values for the support parameter a, we plot in Fig. 4 (left) the transmission coefficients (top) 
and resonances in the lower complex plane (bottom) in corresponding colors. For each potential, 
the location of the resonant frequency, ky t is indicated by a black cross (+) in the transmission 
diagram. The resonances were computed by solving the associated quadratic eigenvalue problem 
using MatScat [3]. 

Remark 6.3. As in Remark 6.1, let q(x) be a periodic potential and qu(x) = <l( x )l\-M,M]( x )- As 
M t oo, the resonances of —d% + qM will converge to the spectrum of —d% + (/oo [2, 17]. 

In Fig. 4 (left), we see from both the transmission coefficient and the resonances that a gap is 
opening in the spectrum as a f oo, supporting Conjecture 1.1, that qoo(x) is a periodic potential 
with a localized defect. 

Remark 6.4. For large support parameter a, a narrow spike forms in the transmission coefficient 
for a value k within the spectral gap of the limiting operator and a resonance lies nearby. In 
the limit that a f oo, this resonance converges to a point eigenvalue within the spectral gap. For 
periodic potentials with a localized defects, such defect eigenvalues exist [16, 9, 10, 27]. Our V op t are 
qualitatively similar to the class studied in [16]. Note that the spike in the transmission coefficient 
in Fig. 4 (left) appears to lie near the resonant frequency, but at a distinct value. 

In Fig. 4 (right) we plot k vs. the transmission coefficient |iy(/c)| 2 for the color-corresponding 
potentials obtained by varying \x (the forcing frequency for which the optimization is performed) 
in Fig. 2 (right). For each value of //, the resonant frequency lies in a spectral gap for each value 
of \x and there appears only to be a single gap. 

Remark 6.5. The Schrodinger operator H q = —d% + q with one specified spectral gap is unique 
and can be explicitly written in terms of Jacobi elliptic functions [15]. These are called one-gap 
potentials. Using the transmission coefficient plots in Fig 4 (right) to estimate the position of the 
spectral gap, we find that the corresponding one-gap potential has period comparable to that of 
the periodic tail of the potentials given in Fig. 2(right). 

This suggests a good heuristic for finding potentials with small values of r[V]: Start with 
a localized potential well supporting a single bound state. Then, create a low density of states 
at k = y/Xy + A* by adding a truncated one-gap potential with appropriate spectral gap. If the 
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Figure 5: With /3 = V as in Eq. (4.11), we plot locally optimal potentials for varying values of 
support a and forcing frequency \x. 



potential added has small amplitude, then this heuristic is equivalent to adding a cosine or Mathieu 
potential with frequency given by the Bragg relation. 

6.5 Optimizing T with (3 = V as in Eq. (4.11) 

Here we study the case where (3 = V as in Remark 1.2, Eq. (4.11), and Corollary 4.13. 

In Fig 5 (left), we take \x = 2 and plot locally optimal potentials for 4 different values of the 
support, a. The values of V are given in the following table. 



a 


4 


8 16 


32 


r 


8 x i(r 13 


3 x 1(T 13 2 x 10 -13 


2 x i(r 12 



In Fig. 5 (right), we take a = 32, and plot locally optimal potentials for 4 values of forcing 
frequency fi. The values of T are given in the following table. 





2 


3 


4 


5 


r 


2 x i(r 12 


1 x 1(T 12 


2 x kt 13 


2 x KT 15 



As noted in Remark 4.1, the solution of the potential design problem is not guaranteed to be 
unique, since the admissible set is non-convex. Regarding Conjecture 1.1 on the character of the 
limit of optimizers, V op t t a as a tends to infinity, since for j3{x) = V(x) and a = oo, the functional 
V I—)- r[V] is invariant under the transformation V(x) \— > V(x + Xq), we could expect convergence 
to V op t t oo(x), a localized perturbation of a periodic potential, only modulo translations. 



6.6 Time dependent simulations 

For a locally optimal potential of the potential design problem, (5.1), we independently verify that 
the potential supports a very long-lived metastable state by conducting time-dependent simulations. 
See Section 5 for a discussion of the numerical method. We set Vmit = —Asech 2 (Bx) for suitably 
chosen A, B and take V op t to be a locally optimal solution to the PDP (5.1) with (3 = 1[_2,2]> i 1 = 2, 
and a = 12 (same parameter choice as in Section 6.1). We then solve the parametrically forced 
Schrodinger Eq. (1.3) with e = 1 until t = 40 with initial conditions given by the ground state of 
Hy for the two potentials, i.e. 4> e (0) = tpv op t an d </> e (0) = ipv ini f I n Figs. 1 (left) and 1 (center) we 
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Figure 6: For the two potentials in Sec. 6.6 and Fig. 1, and cancellation potential in Fig. 3(right), 
we plot the time evolution ^(t, x), governed by Eq. (1.3) with initial condition taken to be 
ipy + noise. The simulation was performed on a spatial domain [—60, 60] with absorbing boundary 
conditions. 

plot V, /3, and ipy for the two potentials. In Figure 1 (right), we plot t vs. \(<p e (t, V'vC")) | 2 > ^ ne 
square modulus of the projection of the wave function onto the bound state for the two potentials. 

6.7 Filtering study 

For the same potentials studied in Sec. 6.6 and Fig. 1 plus the one studied in Fig. 3 (right), we 
conduct the following experiment. We consider the time evolution of Eq. (1.3) until time t = 50 
with initial condition taken to be ipy + noise. The noise is taken to be normally distributed random 
numbers generated using Matlab's randn function for each point in the interval [—a, a]. The initial 
condition is then normalized so that (c/> e (0), ipy) = 1. The results are plotted in Fig. 6. We find 
that for a non-optimized potential, the final state of the system is nearly zero. While for the locally 
optimal potential, the bound state emerges as the final state. In the central panel of Figure 6 we 
see convergence to the projection of the initial condition onto the bound state of Hy A ; see central 
panel of Figure 1. In the right panel of Figure 6 we see convergence to the projection of the initial 
condition onto the bound state of Hy B . 

This study suggests that such a device could be used as a filter to select a particular spatial 
mode profile. For these potentials, the system behaves as a mode-selecting waveguide, preserv- 
ing the discrete components of the initial condition, while radiating the continuous components. 
Alternatively, this study demonstrates the robustness of ipy opt to large fluctuations in the data. 

7 Discussion and conclusions 

Scattering loss, a conservative loss mechanism, is often a limiting factor in the performance of many 
engineered devices. Therefore, there is great interest in finding structures with low scattering loss- 
rate. Loss can occur due to parametric or nonlinear time-dependent perturbations which couple 
an ideally isolated state to an environment. We consider a model of a bound state supported by a 
potential, V, subject to a time-periodic and spatially localized "ionizing" perturbation. The rate 
of scattering loss, r[V], due to coupling of the bound state to radiation modes is given by Fermi's 
Golden Rule, which depends on the potential V. Using gradient-based optimization methods we 
find locally optimal structures with much longer-lived bound states. These potentials appear to be 
truncations of smooth periodic structures with localized defects. This approach can be extended 
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to the wide class of problems presented in the introduction. 



A Computation of gradients / functional derivatives 
A.l Proof of Prop. 4.5b, gradient of r[V] 

Proof. Here we compute the individual terms in Eq. (4.5a) and then assemble below. 
Computation of and s ^ . We compute 

[#] = e v ±)W^v,ev±} + cc. (A.la) 

^-K^(/3Vy,ey±)(/3ey±,^). (A.lb) 



difjy lQky 



Similarly, 



r ' V [Se v ] = ~M(e v ±,pip v )(/3ip v ,6e v ). (A.2) 



6ey± 8ky 



Computation of ^f- and ji^r- Taking variations of Hyijjy = Xyipy we find that 

(Hy — Xy)di/jy = — (SVifrv — SXyipv)- 
Multiplying by ip and integrating, we obtain 5Xy = (ipy,5Vipy), 

w=# (A - 3) 



and 



(Hy — Xv)Sipv = _ (SVipv — (Vv> SVtpy)ipy) 

= -pL Wvsv] . 



where is the orthogonal projection onto the space spanned by tpy. Since Hy supports only a 
operator 



single bound state, P^ y = P c . The solution of this equation can be written in terms of the resolvent 



S -^[SV] = SiPv = -R v (^)P c [ip v 6V] (A.4) 

Computation of Se ^ v ± and ^f- We note that the distorted plane waves ey± can be expressed 

e v± (x,k v ) = e ±lkvx -<p v ±(x,ky) (A.5) 

where (j)y± satisfies the following equation with outgoing boundary conditions 

(Hy - k v )4>y ± = Ve ±lkvx xefl = [-L, L] (A.6a) 
V(f)y± • n = iky(j)y± x 6 dfl. (A. 6b) 
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Taking variations of Eq. (A. 6), we obtain 

(H v - k v )8<j> v ± = SVe v ± + (2k v <j) V± ± ixVe ±lkvX 
V54>v± • n — iky6<f>y± = i5kv4>v±- 



5ky 



(A.7a) 
(A.7b) 



Recalling ky = Xy + \i and Eq. (A. 3) we find that Sky 

5k v _ ip v 



1&W\ = or equivalents 



(A. 



Equation (A. 7) is a forced equation for 8<j)y± with a unique solution since there is no nontrivial, 
outgoing solution to the homogenous equation [34]. The general solution of Eq. (A. 7a) is 



5V 



[5V] = a±ey+ + (3±ey- + Rv(kv) bVey± + ( 2kycj)y± ± ixVe 



Aiikvx 



Sky 



where a, f3 are constants. Matching boundary conditions in Eq. (A. 7b) and recalling that Ry is 
the outgoing resolvent, we obtain 



a± 



5ky 
2ky 

5ky 

2ky~ 



ikyL 



4>V±{-L)e 
4> v ±(L)e^ L 



so that 

5<t>v± , 
SV 



[SV] = ( R v (k v ) \2k v 4> v ± ± ixVe ±lkvl 



ik v L 



2k v 



l )±(-L)e v + + <f)±(L)e V -) 



Sky 
W 



+ R v (k v )[SVe v ±]. 
Now using Eq. (A. 5) we find that 

Se v ± = ±ix5k v e ±lkvx - 5<j>v±- 



[SV] 
(A.9) 



(A.10) 



Computation of Terms in Eq. (4.5a). Using Eqs. (A.l) and (A. 4) we obtain for the first 
term in Eq. (4.5a) 



ST 

Sipy 



[SV] 



(- 



1 



8k v 



ifoR V (y/\ V )P t 



± 



,5V) 



(A.ll) 



where we have used the fact that the operator Ry(\fX)P c : L 2 —¥ L 2 is symmetric. 

The second term of Eq. (4.5a) can be written using Eqs. (A. 2), (A. 8), (A.9), and (A. 10) 



Se v ± 



8ey±(-,ky) 
SV 



[SV] 



1 

1 



M{/3^ v ,e v ±){^ v ,A ± 5k v - R v (k v ) [6Ve v ±]) (A.12a) 



2k v 



(Mv, A ± )i> 2 v - W±Rv{k v ) Wv] , SV) 
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where A± is given in Eq. (4.6) and we have again used the fact that Ry is a symmetric operator. 
Using Eq. (A. 8), the third term of Eq. (4.5a) is given by 

From Eqs. (4.5a), (A. 11), (A. 12), and (A. 13) and the Riesz representation theorem we obtain 
Eq. (4.5b) as desired. □ 

A.2 Proof of Prop. 4.3, gradient of W v (0) 
Proof. Denoting f(x) = ^p-[5V{y)\, we fix x and compute 

W(x) = T]+T]_ + rj + r]'_ — rjifj- — T],r]—. (A. 14) 

To compute rj±, we take variations of Eq. (4.3) to obtain 

Hy5r]± = —5Vrj± 
lim d x 5rj± = 0. 

x— >±oo 

Using the variation of parameters formula, we find 

Sfj-\- (x) f 
V±( x ) = §v [8V] = - q(x,y)5V(y)r]±(y)dy 

where 

q ( x y \ = fv-Wv+iy) x ~y 

w \n+{x)r}-{y) x >v 

is the Green's function. Differentiating we find 

r/±(x) = - d x q(x,y)5V(y)r]±(y)dy. 



We now break Eq. (A. 14) into 2 parts: W = W\ + W2 where W\ = f* o0 *dy, W2 = f£°*dy, and 
the integrand is given by 

* = -6V(y) (q(x,y)r] + (y)ri'_(y) + r] + (x)d x q(x, y)r/_ (y) - d x q(x,y)r] + (y)r]^(x) + Tj' + (x)q(x,y)r]-(y)) 
We then evaluate 

Wx = -~J (rj + (x)ri- (y)y+ {y)rf_ (x) + r? + (x)t)' + (x)t?_ (y)r?_ (y) 

- r]' + (x)r)- {y)r}+ (y)r?_ (y) - rj' + (x)rj + (x)r]-(y)r]-(y) ) 5V(y) dy 

v+(y)v-(y) sv (y) d y 
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where the underlined terms cancel and 

If 00 / 

W 2 = ~— J [ T]-(x)r] + (y)T] + (y)r]'_(x) + rj+ {x)r]_ (x)rj+ (y)rj- (y) 

- yL(x)rj+{y)ri + (y)r]-(x) - ■q l + {x)r}-{x)-q + {y)r]-{y)^8V{y) dy 

poo 

= -/ V+(y)rj-(y)5V(y)dy. 

Thus 

W = - J rj + (y)rt-(y)SV(y)dy 
and the result follows. Note that W(x, y) is constant in x as expected. □ 
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